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Abstract 

The atmospheres of planets (including Earth) and the outer layers of 
stars have often been treated in radiative transfer as plane-parallel me- 
dia, instead of spherical shells, which can lead to inaccuracy, e.g. limb 
darkening. We give an exact solution of the radiative transfer spe- 
cific intensity at all points and directions in a finite spherical medium 
having arbitrary radial spectral distribution of: source (temperature), 
absorption, emission and anisotropic scattering. The power and ef- 
ficiency of the method stems from the spherical numerical gridding 
used to discretize the transfer equations prior to matrix solution: the 
wanted ray and the rays which scatter into it both have the same 
physico-geometric structure. Very good agreement is found with an 
isotropic astrophysical benchmark (Avrett &; Loeser, 1984). We intro- 
duce a specimen arbitrary forward-back-side phase scattering function 
for future comparisons. Our method directly and exactly addresses 
spherical symmetry with anisotropic scattering, and could be used to 
study the Earth's climate, nuclear power (neutron diffusion) and the 
astrophysics of stars and planets. 

PACS numbers: 42.68.Ay 95.30.Jx 97.10.Ex 97.10.Fy 

1 Introduction 

Radiative transfer grew largely from astrophysics [1,2] and later moved to 
terrestrial problems of atmosphere [3], hot gas and solids [4,5] and nuclear 
power [6] . Many of these problems dealt with the outer layers of a spherical 
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medium but, for historical, technical and other reasons, these were treated 
as plane parallel media [1], and the theoretical results were often very satis- 
factory. Nevertheless, there presently appears to be no ab initio direct sys- 
tematic accurate accepted and used method for calculating radiative transfer 
in a a spherical, emitting, absorbing and anisotropically scattering medium. 
We now present such a calculation, but firstly recall the 'general' transfer 
equation (RTE): 

dl(r s) 

^ = -J(r,s)[a(r,s) +ff(r,s)] + a(r, s)5(r, s) 

+^T^ / l'MF(r,s,s')dn(v,s') (1) 

Here I(r, s) is the (specific) intensity of photons inside the medium at 
point r and detector direction s, B(r, s) is the photon source function, a(r, s) 
is the absorption coefficient, cr(r, s) is the scattering coefficient, and a frac- 
tion (part) F(r, s, s')c?J?(r, s') of the non-detector intensity I'(r,s') is scat- 
tered into the detector direction I(r, s); the scattering phase function is 
F(r, s,s') and dQ(r,s') the solid angle around direction s'. 

Before analysing spherical symmetry, we note that (1) has the form of 
a many-particle transport equation: radiative transport for photons [2,3], 
neutron diffusion for nuclear reactions [6] and the Boltzmann equation for 
gas molecules [7]. 



2 Analysis 

We calculate the intensity by numerical solution of the of the sphere gridded 
as in Fig.l. Spherical shells m = 1..M between adjacent radii r m contain 
abritrary constant values of source B m , absorption a m and scattering a m . 
With mid-radii 

r m = \{r m + r m +i), (2) 

we define the intensity 7 m j as a ray tangential to f m with j shells to the left. 
We introduce the standard albedo uj = a / '(a + a) and incremental optical 
depth 

^Tm,j = ^ s m,j(.Ctm,j ~\~ 0~m,j)i (3) 
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where As m j = ab is the physical length step j along I m j. Then the general 
(1) becomes 

. t mm 2m' — 1 

^ = -J m , + (i - c mJ )B mJ + ^ E E l <f C^^r 

u m'=l j'=l 

(4) 

The logic term log m ^ says which of the many I' m , ■, physically scatter at 
and along I m j and, together with step As m j, phase F™,'- 7 and solid-angle 
Afi m / j/, are derived in the Appendix, all as functions of the the shell radii 

We now make an essential observation that in (4) the set of detector 
intensities {I m ,j} is identical to the the set of scattering intensities {I' m , ■,}: 
this is visually obvious from Figure l(a,b), and would not necessarily be 
so had we chosen a different geometric gridding. Thus 1 = 1', and with 
Aim,j = Im,j — Im,j-i, we can simply rearrange (4) as 



AO 



1 + Ar mj 



J m,i — -, , a ~r 



mm 2m' -1 XX i m V,j' p m J 1 » J , , 

^ 1 + Ar mi - m ' ,j ' ( ' 

m'=l j'=l mj 

In (5) we can count the integer index pairs (m,j) and (m',f) by single 
integers i and i' = 1..M 2 : 

i = i(m,j) = (m — l) 2 + j 

i' = i( m ', f) = (rri - l) 2 + j' (6) 

and (5) becomes 

I i = S i + J2 K i,i' I i' (7) 

v 

where the source vector S is 



Bm,j ATjtjJ ( 1 ^m,j 
T m ,j 

and the scattering kernel matrix if is 



*i - ^(mj) " 1 + Ar?? . 



-"-1,1' — -"■i(m,j),i(m'j') — -, , A • l y J 
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The matrix equation (7), I = S + KI, has solution / = (1 — K)~ l S 
given by a standard numerical linear algebra solver, e.g. BLAS. An Jj = 
Itimj) = Im,j is t ne intensity at any point with radius f m _i m _ji and at angle 
sin -1 (f m /f m _| m _j|) to it. Matrix solution / of (7) gives M 2 exact values of 
I m ,j, whether or not we need them all. The I m ,2m-i are surface intensities, 
as detected in vacuo at distance R> n, the others are internal intensities. 

3 Results 

We calculate the astrophysical benchmark of Avrett & Loeser '84 [8] by our 
method. Their results are the basis of other comparisons [9-14]. We have 
implemented our analysis §2 as computer code in both Fortran and Math- 
cad, giving identical numerical results. The results are for the directionally 
averaged intensity J m in shell m, 

M 2m'-l 

J ™ = 4^ E E 5 ^ m ' ~ W -i'\)[l + 5{m',j')]I m , r ^ mj (10) 

m'=l j'=l 

See [8] for details of the code inputs {M, r m ,B m , a m , a m ,F = 1}, isotropic 
scattering. Our results are in (very) good agreement with [8], see Fig. 2, 
which attest to the validity and value of our method. 

4 Discussion and applications 

• All quantities so far are implicictly monochromatic, at photon fre- 
quency v. If shell m! is in local thermodynamic equilibrium (LTE) 
then it may be assigned a temperature T m > , and the source function is 
the Planck black-body function: 

B = = ^ 

mj ^ exp(W^ m _| m _,|)-l 

• The ratio f:s:b determines the forward, side and back strength of the 
anisotropic scattering eglipsoid, see the Appendix. We have calculated 
the intensity 7 m j in an homogenous B = 1 u = 0.5 R = 1 M = 5 
shell sphere for the cases (f,s,b) = (2,2,2), (1,2,3), (1,3,2), (2,1,3), 
(2,3,1), (3,2,1), (3,1,2). These give a quantitative feel for the effects of 
anisotropic scattering, and are available on request from the author. 
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• We can simulate a planetary atmosphere (e.g. Earth) by h shells 
r\..Th as gas and shells r^\..M as solid, with shell/surface having 
emissivity e (and thus reflectivity 1 — e). A planetary atmosphere is 
typically 'thin', r\ — <C vm, and has 'high' absorption Oi so ud S> a gas 
in the solid core. Rays I m j with m < h are 'limb' rays which do 
not intercept the solid planet and are not easily modelled in a plane 
parallel situation, but are rays like any others for our method. If the 
planet receives the Sun's intensity I so i a r from left to right, then we 
have boundary conditions 

4ij — l'rn,j=0 — Irnfi — ^solari (H) 

with day in the left hemisphere and night in the right. 

Realistic code implementation of the the above planetary physics, es- 
pecially our calculation with a total planetary volume mainly of high-a 
solid and its surface, needs careful comparison to present methods to 
judge if the differences between our 'direct-sphere' method and con- 
ventional methods (plane-parallel) are of value and advantage. 



Appendix 

We here give expressions for log™ ^ , As m j, ACl m /ji and F™;-*; all derive 
from elementary spherical/circular geometry. 

Scattering 'selection rules'. From Figure 1(b) we see that the detec- 
tor ray I m j is on shell number 

rrij = m — \m — j\ m j = m2m-j (12) 

and thus only (scattering) rays I rn iji with m'-, = rrij can scatter into it. 
With (10) this means 

log^f =5(m' j „m j )[l + 5(f,m')} (13) 



Physical step length. From Figure 1(c), this is ab = 



As 



\J r m-\m-j\ rfn \J r m-\m-j-l\ ^ if < J < 
2 \/ r m ~ f2 m if j = TH 



rn 



(14) 



As 



m,2m—j 



if 2m > j > m 
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Solid angle. A cone of semi-angle 9 subtends solid angle f2(0) = 2tt(1 — 
cos 6) and, see Figure 1(c) which is a central circular section of the sphere, 
ray I m j is on the cone fl(AjA') whose surface Aj is swept out by revolution 
about axis Oj, and contains the tangent sphere radius f m : the solid angle 
of I m j in (4) is then the conical shell between adjacent cones Q(BjB') and 
£l(CjC), containing tangent spheres r m and r m+ \: 

An m;j = n(BjB>) - n(CjC) 



= 2tt < 



rn+ll m-\m-j\ 



r 2 Jf 2 



m— \m— j\ 



ii 1 



.2 1=2 

rn+ll m-\rn-j\ 



m-j\ 



if m = M and j = m 
if m = M and j / m 

if m < M and j = m 

if m < M and j ^ m 
(15) 



Scattering phase function. The analysis so far directly covers isotropic 
scattering, F = 1, and we now turn to the detail of anisotropic scattering. 
For most scattering media of interest (random, homogenous, not 'oriented' 
or crystalline) the strength of anisotropic scattering depends only on the 
angle 7 = J^^j between the detector ray I Pt m,j an d the scattering ray 
Ip',m',j''i i n Fig- l( c ) with polar axis z = jO and x on the paper surface, 
(j) p = 2ir(p — 1)/P is the azimuthal angle of I Pt m,j on the detector cone m 
and likewise </y is the azimuthal angle of I p ',m',j' on t ne scattering cone m'. 
Then, with sin0 m;j = r m /f m _\ m _j\, 

cos ^p,m,'j = cos °m,j cos 6» m /j + sin m ,j sin m / ;j cos(0 p - p ,) (16) 
Since symmetry says the magnitude of intensity I P , m ,j is independent of 4> p , 
we choose <p p = and can thus write ip,m,j = ^m,j and cos 7^ ^ ■ = cos 7^'™ - 

The scattering strength is given by some function F'{p/) defining a non- 
spherical surface. For illustration we choose surface F' to be egg-shaped, 
being the smooth union of two different half-ellipsoids of revolution, an eglip- 
soid, defined by forward /, side s and back b major/minor axes: 

f'(v) = J N v> b > s W v2 f~ 2 + g EZEf if 1 > ^ > , m 

{) \ N(f,b,s)/^v 2 b- 2 + (l-v 2 )s- 2 if-l<v<0 1 ' 
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where N ensures phase normalisation, J F'dVt = Att. Thus a fraction F'^^'J 1 ) 

of I p \ m 'j' is scattered along I m j (=Ip, m ,j), and phase function F™; J is the 
sum of all p' = 1, 2, ..P such contributions: 

P '=i 
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Figure 1: Geometric spherical gridding of our analysis, §2; (a) and (b) show 
that the set {/} of detector rays is identical to the scattering set {/'}, and 
(c) shows constructions for step As m j, solid angle AQ m j and polar axes 
xyz for the azimuthal <p p of anisotropic scattering. 



Figure 2: Comparison of benchmark [8] (circles), to our results (lines) §3. 
The M =14 shell sphere has optical density (depth) r oc 1/r. 



